{
    "cells": [
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "# Predicting hourly wind speed with JFK airport weather data"
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "if not ('sc' in locals() or 'sc' in globals()):\n    print('It seems you are not running in a IBM Watson Studio Apache Spark Notebook. You might be running in a IBM Watson Studio Default Runtime or outside IBM Waston Studio. Therefore installing local Apache Spark environment for you. Please do not use in Production')\n    \n    from pip import main\n    main(['install', 'pyspark==2.4.5'])\n    \n    from pyspark import SparkContext, SparkConf\n    from pyspark.sql import SparkSession\n\n    sc = SparkContext.getOrCreate(SparkConf().setMaster(\"local[*]\"))\n    \n    spark = SparkSession \\\n        .builder \\\n        .getOrCreate()",
            "execution_count": null,
            "outputs": []
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Let\u2019s start by downloading the dataset and creating a dataframe. This dataset can be found on DAX, the IBM Data Asset Exchange and can be downloaded for free.\n\nhttps://developer.ibm.com/exchanges/data/all/jfk-weather-data/"
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "# delete files from previous runs\n!rm -f jfk_weather*\n\n# download the file containing the data in CSV format\n!wget http://max-training-data.s3-api.us-geo.objectstorage.softlayer.net/noaa-weather/jfk_weather.tar.gz\n\n# extract the data\n!tar xvfz jfk_weather.tar.gz\n    \n# create a dataframe out of it by using the first row as field names and trying to infer a schema based on contents\ndf = spark.read.option(\"header\", \"true\").option(\"inferSchema\",\"true\").csv('jfk_weather.csv')\n\n# register a corresponding query table\ndf.createOrReplaceTempView('df')",
            "execution_count": 1,
            "outputs": [
                {
                    "output_type": "stream",
                    "text": "Waiting for a Spark session to start...\nSpark Initialization Done! ApplicationId = app-20200407032048-0000\nKERNEL_ID = 0d4d9307-8faa-4f9d-9c7a-b7f885034dab\n--2020-04-07 03:20:52--  http://max-training-data.s3-api.us-geo.objectstorage.softlayer.net/noaa-weather/jfk_weather.tar.gz\nResolving max-training-data.s3-api.us-geo.objectstorage.softlayer.net (max-training-data.s3-api.us-geo.objectstorage.softlayer.net)... 67.228.254.196\nConnecting to max-training-data.s3-api.us-geo.objectstorage.softlayer.net (max-training-data.s3-api.us-geo.objectstorage.softlayer.net)|67.228.254.196|:80... connected.\nHTTP request sent, awaiting response... 200 OK\nLength: 2575759 (2.5M) [application/x-tar]\nSaving to: 'jfk_weather.tar.gz'\n\n100%[======================================>] 2,575,759   --.-K/s   in 0.03s   \n\n2020-04-07 03:20:52 (80.8 MB/s) - 'jfk_weather.tar.gz' saved [2575759/2575759]\n\n./._jfk_weather.csv\njfk_weather.csv\n",
                    "name": "stdout"
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "The dataset contains some null values, therefore schema inference didn\u2019t work properly for all columns, in addition, a column contained trailing characters, so we need to clean up the data set first. This is a normal task in any data science project since your data is never clean, don\u2019t worry if you don\u2019t understand all code, you won\u2019t be asked about it. "
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "import random\nrandom.seed(42)\n\nfrom pyspark.sql.functions import translate, col\n\ndf_cleaned = df \\\n    .withColumn(\"HOURLYWindSpeed\", df.HOURLYWindSpeed.cast('double')) \\\n    .withColumn(\"HOURLYWindDirection\", df.HOURLYWindDirection.cast('double')) \\\n    .withColumn(\"HOURLYStationPressure\", translate(col(\"HOURLYStationPressure\"), \"s,\", \"\")) \\\n    .withColumn(\"HOURLYPrecip\", translate(col(\"HOURLYPrecip\"), \"s,\", \"\")) \\\n    .withColumn(\"HOURLYRelativeHumidity\", translate(col(\"HOURLYRelativeHumidity\"), \"*\", \"\")) \\\n    .withColumn(\"HOURLYDRYBULBTEMPC\", translate(col(\"HOURLYDRYBULBTEMPC\"), \"*\", \"\")) \\\n\ndf_cleaned =   df_cleaned \\\n                    .withColumn(\"HOURLYStationPressure\", df_cleaned.HOURLYStationPressure.cast('double')) \\\n                    .withColumn(\"HOURLYPrecip\", df_cleaned.HOURLYPrecip.cast('double')) \\\n                    .withColumn(\"HOURLYRelativeHumidity\", df_cleaned.HOURLYRelativeHumidity.cast('double')) \\\n                    .withColumn(\"HOURLYDRYBULBTEMPC\", df_cleaned.HOURLYDRYBULBTEMPC.cast('double')) \\\n\ndf_filtered = df_cleaned.filter(\"\"\"\n    HOURLYWindSpeed <> 0\n    and HOURLYWindSpeed IS NOT NULL\n    and HOURLYWindDirection IS NOT NULL\n    and HOURLYStationPressure IS NOT NULL\n    and HOURLYPressureTendency IS NOT NULL\n    and HOURLYPrecip IS NOT NULL\n    and HOURLYRelativeHumidity IS NOT NULL\n    and HOURLYDRYBULBTEMPC IS NOT NULL\n\"\"\")",
            "execution_count": 2,
            "outputs": []
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "We want to predict the value of one column based of some others. It is sometimes helpful to print a correlation matrix. "
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "from pyspark.ml.feature import VectorAssembler\nvectorAssembler = VectorAssembler(inputCols=[\"HOURLYWindSpeed\",\"HOURLYWindDirection\",\"HOURLYStationPressure\"],\n                                  outputCol=\"features\")\ndf_pipeline = vectorAssembler.transform(df_filtered)\nfrom pyspark.ml.stat import Correlation\nCorrelation.corr(df_pipeline,\"features\").head()[0].toArray()",
            "execution_count": 3,
            "outputs": [
                {
                    "output_type": "execute_result",
                    "execution_count": 3,
                    "data": {
                        "text/plain": "array([[ 1.        ,  0.25478863, -0.26171147],\n       [ 0.25478863,  1.        , -0.13377466],\n       [-0.26171147, -0.13377466,  1.        ]])"
                    },
                    "metadata": {}
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "vectorAssembler2 = VectorAssembler(inputCols=[\"HOURLYWindSpeed\",\"HOURLYPressureTendency\"],\n                                  outputCol=\"features\")\ndf_pipeline2 = vectorAssembler2.transform(df_filtered)\nCorrelation.corr(df_pipeline2,\"features\").head()[0].toArray()",
            "execution_count": 10,
            "outputs": [
                {
                    "output_type": "execute_result",
                    "execution_count": 10,
                    "data": {
                        "text/plain": "array([[ 1.        , -0.01324305],\n       [-0.01324305,  1.        ]])"
                    },
                    "metadata": {}
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "As we can see, HOURLYWindSpeed and HOURLYWindDirection correlate with 0.06306013 whereas HOURLYWindSpeed  and HOURLYStationPressure correlate with -0.4204518, this is a good sign if we want to predict HOURLYWindSpeed from HOURLYWindDirection and HOURLYStationPressure.\nSince this is supervised learning, let\u2019s split our data into train (80%) and test (20%) set."
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "splits = df_filtered.randomSplit([0.8, 0.2])\ndf_train = splits[0]\ndf_test = splits[1]",
            "execution_count": 4,
            "outputs": []
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Again, we can re-use our feature engineering pipeline"
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "from pyspark.ml.feature import StringIndexer, OneHotEncoder\nfrom pyspark.ml.linalg import Vectors\nfrom pyspark.ml.feature import VectorAssembler\nfrom pyspark.ml.feature import Normalizer\nfrom pyspark.ml import Pipeline\n\nvectorAssembler = VectorAssembler(inputCols=[\n                                    \"HOURLYWindDirection\",\n                                    \"ELEVATION\",\n                                    \"HOURLYStationPressure\"],\n                                  outputCol=\"features\")\n\nnormalizer = Normalizer(inputCol=\"features\", outputCol=\"features_norm\", p=1.0)",
            "execution_count": 5,
            "outputs": []
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Now we define a function for evaluating our regression prediction performance. We\u2019re using RMSE (Root Mean Squared Error) here , the smaller the better\u2026\n\n"
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "def regression_metrics(prediction):\n    from pyspark.ml.evaluation import RegressionEvaluator\n    evaluator = RegressionEvaluator(\n    labelCol=\"HOURLYWindSpeed\", predictionCol=\"prediction\", metricName=\"rmse\")\n    rmse = evaluator.evaluate(prediction)\n    print(\"RMSE on test data = %g\" % rmse)",
            "execution_count": 6,
            "outputs": []
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Let\u2019s run a linear regression model first for building a baseline.\n\n"
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "#LR1\n\nfrom pyspark.ml.regression import LinearRegression\n\n\nlr = LinearRegression(labelCol=\"HOURLYWindSpeed\", featuresCol='features_norm', maxIter=100, regParam=0.0, elasticNetParam=0.0)\npipeline = Pipeline(stages=[vectorAssembler, normalizer,lr])\nmodel = pipeline.fit(df_train)\nprediction = model.transform(df_test)\nregression_metrics(prediction)",
            "execution_count": 9,
            "outputs": [
                {
                    "output_type": "stream",
                    "text": "RMSE on test data = 5.53835\n",
                    "name": "stdout"
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "#GBT1\n\nfrom pyspark.ml.regression import GBTRegressor\ngbt = GBTRegressor(labelCol=\"HOURLYWindSpeed\", maxIter=100)\npipeline = Pipeline(stages=[vectorAssembler, normalizer,gbt])\nmodel = pipeline.fit(df_train)\nprediction = model.transform(df_test)\nregression_metrics(prediction)",
            "execution_count": 8,
            "outputs": [
                {
                    "output_type": "stream",
                    "text": "RMSE on test data = 5.10736\n",
                    "name": "stdout"
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Now let\u2019s switch gears. Previously, we tried to predict HOURLYWindSpeed, but now we predict HOURLYWindDirection. In order to turn this into a classification problem we discretize the value using the Bucketizer. The new feature is called HOURLYWindDirectionBucketized."
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "from pyspark.ml.feature import Bucketizer, OneHotEncoder\nbucketizer = Bucketizer(splits=[ 0, 180, float('Inf') ],inputCol=\"HOURLYWindDirection\", outputCol=\"HOURLYWindDirectionBucketized\")\nencoder = OneHotEncoder(inputCol=\"HOURLYWindDirectionBucketized\", outputCol=\"HOURLYWindDirectionOHE\")\n",
            "execution_count": 11,
            "outputs": []
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Again, we define a function in order to assess how we perform. Here we just use the accuracy measure which gives us the fraction of correctly classified examples. Again, 0 is bad, 1 is good."
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "def classification_metrics(prediction):\n    from pyspark.ml.evaluation import MulticlassClassificationEvaluator\n    mcEval = MulticlassClassificationEvaluator().setMetricName(\"accuracy\") .setPredictionCol(\"prediction\").setLabelCol(\"HOURLYWindDirectionBucketized\")\n    accuracy = mcEval.evaluate(prediction)\n    print(\"Accuracy on test data = %g\" % accuracy)",
            "execution_count": 12,
            "outputs": []
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Again, for baselining we use LogisticRegression."
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "#LGReg1\n\nfrom pyspark.ml.classification import LogisticRegression\nlr = LogisticRegression(labelCol=\"HOURLYWindDirectionBucketized\", maxIter=10)\n#,\"ELEVATION\",\"HOURLYStationPressure\",\"HOURLYPressureTendency\",\"HOURLYPrecip\"\n\nvectorAssembler = VectorAssembler(inputCols=[\"HOURLYWindSpeed\",\"HOURLYDRYBULBTEMPC\"],\n                                  outputCol=\"features\")\n\npipeline = Pipeline(stages=[bucketizer,vectorAssembler,normalizer,lr])\nmodel = pipeline.fit(df_train)\nprediction = model.transform(df_test)\nclassification_metrics(prediction)",
            "execution_count": 13,
            "outputs": [
                {
                    "output_type": "stream",
                    "text": "Accuracy on test data = 0.692922\n",
                    "name": "stdout"
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "markdown",
            "source": "Let\u2019s try some other Algorithms and see if model performance increases. It\u2019s also important to tweak other parameters like parameters of individual algorithms (e.g. number of trees for RandomForest) or parameters in the feature engineering pipeline, e.g. train/test split ratio, normalization, bucketing, \u2026"
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "#RF1\n\nfrom pyspark.ml.classification import RandomForestClassifier\nrf = RandomForestClassifier(labelCol=\"HOURLYWindDirectionBucketized\", numTrees=30)\n\nvectorAssembler = VectorAssembler(inputCols=[\"HOURLYWindSpeed\",\"HOURLYDRYBULBTEMPC\",\"ELEVATION\",\"HOURLYStationPressure\",\"HOURLYPressureTendency\",\"HOURLYPrecip\"],\n                                  outputCol=\"features\")\n\npipeline = Pipeline(stages=[bucketizer,vectorAssembler,normalizer,rf])\nmodel = pipeline.fit(df_train)\nprediction = model.transform(df_test)\nclassification_metrics(prediction)",
            "execution_count": 14,
            "outputs": [
                {
                    "output_type": "stream",
                    "text": "Accuracy on test data = 0.71895\n",
                    "name": "stdout"
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "#GBT2\n\nfrom pyspark.ml.classification import GBTClassifier\ngbt = GBTClassifier(labelCol=\"HOURLYWindDirectionBucketized\", maxIter=100)\n\nvectorAssembler = VectorAssembler(inputCols=[\"HOURLYWindSpeed\",\"HOURLYDRYBULBTEMPC\",\"ELEVATION\",\"HOURLYStationPressure\",\"HOURLYPressureTendency\",\"HOURLYPrecip\"],\n                                  outputCol=\"features\")\n\npipeline = Pipeline(stages=[bucketizer,vectorAssembler,normalizer,gbt])\nmodel = pipeline.fit(df_train)\nprediction = model.transform(df_test)\nclassification_metrics(prediction)",
            "execution_count": 15,
            "outputs": [
                {
                    "output_type": "stream",
                    "text": "Accuracy on test data = 0.73242\n",
                    "name": "stdout"
                }
            ]
        },
        {
            "metadata": {},
            "cell_type": "code",
            "source": "# rf2\n\nrf2 = RandomForestClassifier(labelCol=\"HOURLYWindDirectionBucketized\", numTrees=10)\n\nvectorAssembler = VectorAssembler(inputCols=[\"HOURLYWindSpeed\",\"HOURLYDRYBULBTEMPC\",\"ELEVATION\",\"HOURLYStationPressure\",\"HOURLYPressureTendency\",\"HOURLYPrecip\"],\n                                  outputCol=\"features\")\n\npipeline2 = Pipeline(stages=[bucketizer,vectorAssembler,normalizer,rf2])\nmodel2 = pipeline2.fit(df_train)\nprediction2 = model2.transform(df_test)\nclassification_metrics(prediction2)",
            "execution_count": 16,
            "outputs": [
                {
                    "output_type": "stream",
                    "text": "Accuracy on test data = 0.721461\n",
                    "name": "stdout"
                }
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "name": "python36",
            "display_name": "Python 3.6 with Spark",
            "language": "python3"
        },
        "language_info": {
            "mimetype": "text/x-python",
            "nbconvert_exporter": "python",
            "name": "python",
            "pygments_lexer": "ipython3",
            "version": "3.6.8",
            "file_extension": ".py",
            "codemirror_mode": {
                "version": 3,
                "name": "ipython"
            }
        }
    },
    "nbformat": 4,
    "nbformat_minor": 1
}